###### R 4.4.3


############# COMBINED NET EFFECT MODEL ########################
reds = reds %>%
  group_by(stateid, districtid, tehsilid, blockid, villageid) %>%
  mutate(vwrXvst_gp1 = ifelse(vwreserv==1 & vstreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwrXvst_gp2 = ifelse(vwreserv==1 & vstreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwrXvst_gp3 = ifelse(vwreserv==1 & vstreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwrXvst_gp12 = ifelse(vwrXvst_gp1==1 | vwrXvst_gp2==1, 1, 0)) %>%
  mutate(vwrXvst_any = ifelse(vwrXvst_gp1==1 | vwrXvst_gp2==1 | vwrXvst_gp3==1, 1, 0)) %>%
  mutate(vwreserv_any = ifelse(vwreserv_gp1==1 | vwreserv_gp2==1 | vwreserv_gp3==1, 1, 0)) %>%
  mutate(vstreserv_any = ifelse(vstreserv_gp1==1 | vstreserv_gp2==1 | vstreserv_gp3==1, 1, 0))
reds_p3_rand = reds %>% filter(nr_short!=1 & period_new==3)


dn.dta <- dn.dta %>%
  mutate(reserved_scst = ifelse(reserved_st==1|reserved_sc==1, 1, 0))



# INTER-CASTE MARRIAGE
controls_dn = c("age + sc + st + fem")
base_dn = paste("marry_other_caste_bin ~ vwreserv*vstreserv")
dn.lm1 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata")

dn.lm1.pval = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lh_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata",
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv")
dn.lm1.pval = dn.lm1.pval$lh$p.value

dn.lm1.pval.st = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lh_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata",
            linear_hypothesis = "vstreserv = vstreserv + vwreserv:vstreserv")
dn.lm1.pval.st = dn.lm1.pval.st$lh$p.value





# INTER-CASTE MARRIAGE - SC/ST
base_dn = paste("marry_other_caste_bin ~ vwreserv*vscstreserv")
dn.lm2 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata")


dn.lm2.pval = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lh_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata",
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")
dn.lm2.pval = dn.lm2.pval$lh$p.value


dn.lm2.pval.stsc = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lh_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata",
            linear_hypothesis = "vscstreserv = vscstreserv + vwreserv:vscstreserv")
dn.lm2.pval.stsc = dn.lm2.pval.stsc$lh$p.value

# EASE OF INTER-CASTE INTERACTIONS
base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"
reds.lm1 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm1.pval = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lh_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata",
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv") 
reds.lm1.pval = reds.lm1.pval$lh$p.value


reds.lm1.pval.st = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lh_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata",
            linear_hypothesis = "vstreserv = vstreserv + vwreserv:vstreserv") 
reds.lm1.pval.st = reds.lm1.pval.st$lh$p.value







# EASE OF INTER-CASTE INTERACTIONS - SC/ST 
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vscstreserv_gp12"
base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"
reds.lm2 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")


reds.lm2.pval = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lh_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")
reds.lm2.pval = reds.lm2.pval$lh$p.value


reds.lm2.pval.scst = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lh_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vscstreserv = vscstreserv + vwreserv:vscstreserv")
reds.lm2.pval.scst = reds.lm2.pval.scst$lh$p.value


# CASTE-BASED CONFLICT
base = "caste_based_inter_family_disputes ~ vwreserv*vstreserv"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_prior + vstreserv_prior + prop_STnow"
reds.lm3 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.st.matched, fixed_effects = ~period_new + districtid, clusters = villageid, se_type = "stata")

reds.lm3.pval = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.st.matched, fixed_effects = ~period_new + districtid, 
            clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv")
reds.lm3.pval = reds.lm3.pval$lh$p.value


reds.lm3.pval.st = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.st.matched, fixed_effects = ~period_new + districtid, 
            clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vstreserv = vstreserv + vwreserv:vstreserv")
reds.lm3.pval.st = reds.lm3.pval.st$lh$p.value









# CASTE-BASED CONFLICT
base = "caste_based_inter_family_disputes ~ vwreserv*vscstreserv"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_prior + vscstreserv_prior + prop_STnow + prop_SCnow"
reds.lm4 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.scst.matched, fixed_effects = ~period_new + districtid, clusters = villageid, se_type = "stata")

reds.lm4.pval = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.scst.matched, fixed_effects = ~period_new + districtid,
            clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")
reds.lm4.pval = reds.lm4.pval$lh$p.value


reds.lm4.pval.scst = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.scst.matched, fixed_effects = ~period_new + districtid,
            clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vscstreserv = vscstreserv + vwreserv:vscstreserv")
reds.lm4.pval.scst = reds.lm4.pval.scst$lh$p.value


# Create dataset of models to loop over
effect_list = list(
  list(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st), outcome = "Marriage", sample = "DN - Entire", model = dn.lm1, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = dn.lm1.pval, coef_dif_p_st = dn.lm1.pval.st),
  list(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_scst), outcome = "Marriage", sample = "DN - Entire", model = dn.lm2, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = dn.lm2.pval, coef_dif_p_st = dn.lm2.pval.stsc),
  list(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1), outcome = "Interaction", sample = "REDS - Entire", model = reds.lm1, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm1.pval, coef_dif_p_st = reds.lm1.pval.st),
  list(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1), outcome = "Interaction", sample = "REDS - Entire", model = reds.lm2, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm2.pval, coef_dif_p_st = reds.lm2.pval.scst),
  list(data = reds.st.matched, outcome = "Conflict", sample = "REDS - Entire", model = reds.lm3, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm3.pval, coef_dif_p_st = reds.lm3.pval.st),
  list(data = reds.scst.matched, outcome = "Conflict", sample = "REDS - Entire", model = reds.lm4, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm4.pval, coef_dif_p_st = reds.lm4.pval.scst))

# Create tables and plots of net effects
plots = lapply(effect_list, FUN = function(x){
  effect_size_plot(data = x$data, outcome = x$outcome, sample = x$sample, model = x$model, 
                   vstvar = x$vstvar, vwvar = x$vwvar, coef_dif_p_w = x$coef_dif_p_w, 
                   coef_dif_p_st = x$coef_dif_p_st)
})

## Net effect tables
effect_tables = bind_rows(plots[[1]] %>% mutate(outcome="Marriage") %>% mutate(quota = "ST"),
                          plots[[2]] %>% mutate(outcome="Marriage") %>% mutate(quota = "ST/SC"),
                          plots[[3]] %>% mutate(outcome="Interaction") %>% mutate(quota = "ST"),
                          plots[[4]] %>% mutate(outcome="Interaction") %>% mutate(quota = "ST/SC"),
                          plots[[5]] %>% mutate(outcome="Conflict") %>% mutate(quota = "ST"),
                          plots[[6]] %>% mutate(outcome="Conflict") %>% mutate(quota = "ST/SC"))
effect_tables = effect_tables %>% 
  mutate(coef_dif_p_chr = ifelse(!is.na(coef_dif_p), paste0("Coef. test p-val=", round(coef_dif_p, 3)), NA_character_))

## Combined plots for REDS sample
m1 = tidy(dn.lm1) %>% mutate(outcome="Marriage")
m2 = tidy(reds.lm1) %>% mutate(outcome="Interaction")
m3 = tidy(reds.lm2) %>% mutate(outcome="Caste Conflict")







####### Interaction effect ######
effect_tables %>% 
  filter(outcome=="Interaction") %>%
  mutate(lab = as.character(lab_order)) %>%
  filter(lab=="Women's Res."|lab=="Combined Women Effect") %>%
  mutate(lab = ifelse(lab=="Combined Women Effect", "Women's Quota \n When Caste Quota=1", lab)) %>%
  mutate(lab = ifelse(lab=="Women's Res.", "Women's Quota \n When Caste Quota=0", lab)) %>%
  mutate(labels = ifelse(lab_order=="Women's Res.", 2, 1)) %>%
  mutate(lab_fac = factor(labels, labels = c("Two-Dimensional\nWomen's Quota\n(Caste Quota=1)", "One-Dimensional\nWomen's Quota\n(Caste Quota=0)"))) %>%
  ggplot(., aes(x = lab_fac, color = quota))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray48")+
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate + 1.96*std.error), size = 1.5,  position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error), size = 3, position = position_dodge(width = 0.4))+
  geom_pointrange(size = 1.5, 
                  aes(y = estimate, shape = quota, 
                      ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error),
                  position = position_dodge(width = 0.4), fill = "white") + 
  geom_text(aes(y = estimate, label = round(estimate,2)), vjust = 2, size = 6,
            position = position_dodge(width = 0.4), show.legend = F)+
  geom_text(aes(y = 0.06, label = coef_dif_p_chr), size = 6, hjust = 0,
            position = position_dodge(width = 0.4), show.legend = F)+
  scale_color_manual(values = c("#e66101", "#b2abd2"), name = "Caste Quota Analyzed")+
  scale_shape_manual(values = c(21, 24), name = "Caste Quota Analyzed")+
  coord_flip()+
  facet_grid(.~outcome, scales = "free")+
  labs(x = "", y = "Estimate (95% CI)")+
  theme_pubr()+
  theme(strip.text = element_text(size = 20), 
        axis.text = element_text(size = 20), 
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 20))
ggsave(paste0(fig.out, "Figure1a.pdf"), 
       width = 10, height = 6.5)

###### Caste quota for interactions ######
effect_tables %>% 
  filter(outcome=="Interaction") %>%
  mutate(lab = as.character(lab_order)) %>%
  filter(lab=="ST Res."|lab=="Combined ST Effect") %>%
  mutate(lab = ifelse(lab=="Combined ST Effect", "Caste Quota\n When Women's Quota=1", lab)) %>%
  mutate(lab = ifelse(lab=="ST Res.", "Caste Quota\n When Women's Quota=0", lab)) %>%
  mutate(labels = ifelse(lab_order=="ST Res.", 2, 1)) %>%
  mutate(lab_fac = factor(labels, labels = c("Two-Dimensional\nCaste Quota\n(Women's Quota=1)", "One-Dimensional\nCaste Quota\n(Women's Quota=0)"))) %>%
  ggplot(., aes(x = lab_fac, color = quota))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray48")+
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate + 1.96*std.error), size = 1.5,  position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error), size = 3, position = position_dodge(width = 0.4))+
  geom_pointrange(size = 1.5, aes(y = estimate, shape = quota, ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error),
                  position = position_dodge(width = 0.4), fill = "white") + 
  geom_text(aes(label = round(estimate,2), y = estimate), vjust = 2, size = 6,
            position = position_dodge(width = 0.4), show.legend = F)+
  geom_text(aes(y = 0.006, label = coef_dif_p_chr), size = 5, hjust = 0,
            position = position_dodge(width = 0.4), show.legend = F)+
  scale_color_manual(values = c("#e66101", "#b2abd2"), name = "Caste Quota Analyzed")+
  scale_shape_manual(values = c(21, 24), name = "Caste Quota Analyzed")+
  coord_flip()+
  facet_grid(.~outcome, scales = "free")+
  labs(x = "", y = "Estimate (95% CI)")+
  theme_pubr()+
  theme(strip.text = element_text(size = 20), 
        axis.text.y = element_text(size = 20),
        axis.text.x = element_text(size = 20),
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 20))

ggsave(paste0(fig.out, "Figure2.pdf"),
       width = 11.6, height = 8.56)


####### Conflict & Marriage effect ######
effect_tables %>% 
  filter(outcome=="Marriage"|outcome=="Conflict") %>%
  mutate(out_fac = ifelse(outcome=="Marriage", 1, 2)) %>%
  mutate(out_fac = factor(out_fac, labels = c("Marriage", "Conflict"))) %>%
  mutate(lab = as.character(lab_order)) %>%
  filter(lab=="Women's Res."|lab=="Combined Women Effect") %>%
  mutate(lab = ifelse(lab=="Combined Women Effect", "Women's Quotas\n With Caste Quotas", lab)) %>%
  mutate(lab = ifelse(lab=="Women's Res.", "Women's Quotas\n Without Caste Quotas", lab)) %>%
  mutate(labels = ifelse(lab_order=="Women's Res.", 2, 1)) %>%
  mutate(lab_fac = factor(labels, labels = c("Two-Dimensional\nWomen's Quota\n(Caste Quota=1)", "One-Dimensional\nWomen's Quota\n(Caste Quota=0)"))) %>%
  ggplot(., aes(x = lab_fac, y = estimate, color = quota))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray48")+
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate + 1.96*std.error), size = 1.5,  position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error), size = 3, position = position_dodge(width = 0.4))+
  geom_pointrange(size = 1.5, aes(y = estimate, shape = quota, ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error),
                  position = position_dodge(width = 0.4), fill = "white") + 
  geom_text(aes(label = round(estimate,2), y = estimate), vjust = 2, size = 6,
            position = position_dodge(width = 0.4), show.legend = F)+
  geom_text(aes(label = coef_dif_p_chr, y = 0.12), size = 4, hjust = 0,
            position = position_dodge(width = 0.4), show.legend = F)+
  scale_color_manual(values = c("#e66101", "#b2abd2"), name = "Caste Quota Analyzed")+
  scale_shape_manual(values = c(21, 24), name = "Caste Quota Analyzed")+
  coord_flip()+
  lims(y = c(-0.15, 0.40))+
  facet_grid(.~out_fac, scales = "free")+
  labs(x = "", y = "Estimate (95% CI)")+
  theme_pubr()+
  theme(strip.text = element_text(size = 20), 
        axis.text.y = element_text(size = 20), 
        axis.text.x = element_text(size = 16),
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 20))
ggsave(paste0(fig.out, "Figure3.pdf"),
       width = 12, height = 7)

####### Combined caste quota ######
effect_tables %>% 
  filter(outcome=="Marriage"|outcome=="Conflict") %>%
  mutate(out_fac = ifelse(outcome=="Marriage", 1, 2)) %>%
  mutate(out_fac = factor(out_fac, labels = c("Marriage", "Conflict"))) %>%
  mutate(lab = as.character(lab_order)) %>%
  filter(lab=="ST Res."|lab=="Combined ST Effect") %>%
  mutate(lab = ifelse(lab=="Combined ST Effect", "Caste Quota\n When Women's Quota=1", lab)) %>%
  mutate(lab = ifelse(lab=="ST Res.", "Caste Quota\n When Women's Quota=0", lab)) %>%
  mutate(labels = ifelse(lab_order=="ST Res.", 2, 1)) %>%
  mutate(lab_fac = factor(labels, labels = c("Two-Dimensional\nCaste Quota\n(Women's Quota=1)", "One-Dimensional\nCaste Quota\n(Women's Quota=0)"))) %>%
  ggplot(., aes(x = lab_fac, color = quota))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray48")+
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate + 1.96*std.error), size = 1.5,  position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error), size = 3, position = position_dodge(width = 0.4))+
  geom_pointrange(size = 1.5, aes(y = estimate, shape = quota, ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error),
                  position = position_dodge(width = 0.4), fill = "white") + 
  geom_text(aes(label = round(estimate,2), y = estimate), vjust = 2, size = 6,
            position = position_dodge(width = 0.4), show.legend = F)+
  scale_color_manual(values = c("#e66101", "#b2abd2"), name = "Caste Quota Analyzed")+
  scale_shape_manual(values = c(21, 24), name = "Caste Quota Analyzed")+
  coord_flip()+
  facet_grid(.~out_fac, scales = "free")+
  labs(x = "", y = "Estimate (95% CI)")+
  theme_pubr()+
  theme(strip.text = element_text(size = 20), 
        axis.text.y = element_text(size = 20),
        axis.text.x = element_text(size = 16),
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 20))
ggsave(paste0(fig.out, "FigureE7.pdf"), 
       width = 12, height = 7)

####### Combined quota (women subsample) ######
# INTER-CASTE MARRIAGE
controls_dn = "age + sc + st"
base_dn = paste("marry_other_caste_bin ~ vwreserv*vstreserv")
dn.lm1 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st) %>%
              filter(fem==1),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata")

dn.lm1.pval = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lh_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st) %>%
              filter(fem==1),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv")
dn.lm1.pval = dn.lm1.pval$lh$p.value


# INTER-CASTE MARRIAGE - SC/ST
base_dn = paste("marry_other_caste_bin ~ vwreserv*vscstreserv")
dn.lm2 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst) %>%
              filter(fem==1),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata")

dn.lm2.pval = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lh_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst) %>%
              filter(fem==1),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata",
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")

dn.lm2.pval = dn.lm2.pval$lh$p.value

# EASE OF INTER-CASTE INTERACTIONS
base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"
reds.lm1 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1) %>%
              filter(fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm1.pval = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lh_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1) %>%
              filter(fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv") 

reds.lm1.pval = reds.lm1.pval$lh$p.value

# EASE OF INTER-CASTE INTERACTIONS - SC/ST 
controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_gp12 + vscstreserv_gp12"
base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"

reds.lm2 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1) %>%
              filter(fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")

reds.lm2.pval = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lh_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1) %>%
              filter(fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")
reds.lm2.pval = reds.lm2.pval$lh$p.value


# # CASTE-BASED CONFLICT
base = "caste_based_inter_family_disputes ~ vwreserv*vstreserv"
controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_prior + vstreserv_prior + prop_STnow"
reds.lm3 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(fem==1), fixed_effects = ~period_new + districtid, clusters = villageid, se_type = "stata")

reds.lm3.pval = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.st.matched %>% filter(fem==1),
            fixed_effects = ~period_new + districtid, clusters = villageid,
            se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv")
reds.lm3.pval = reds.lm3.pval$lh$p.value



# # CASTE-BASED CONFLICT
base = "caste_based_inter_family_disputes ~ vwreserv*vscstreserv"
controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_prior + vscstreserv_prior + prop_STnow + prop_SCnow"
reds.lm4 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(fem==1), fixed_effects = ~period_new + districtid, clusters = villageid, se_type = "stata")


reds.lm4.pval = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.scst.matched %>% filter(fem==1),
            fixed_effects = ~period_new + districtid, clusters = villageid,
            se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")
reds.lm4.pval = reds.lm4.pval$lh$p.value

# Create dataset of models to loop over
effect_list = list(
  list(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st) %>% filter(fem==1), outcome = "Marriage", sample = "DN - Entire", model = dn.lm1, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = dn.lm1.pval, coef_dif_p_st = dn.lm1.pval),
  list(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_scst) %>% filter(fem==1), outcome = "Marriage", sample = "DN - Entire", model = dn.lm2, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = dn.lm2.pval, coef_dif_p_st = dn.lm2.pval),
  list(data = reds.st.matched %>% filter(fem==1) %>% filter(period_new==3 & nr_short!=1), outcome = "Interaction", sample = "REDS - Entire", model = reds.lm1, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm1.pval, coef_dif_p_st = reds.lm1.pval),
  list(data = reds.scst.matched %>% filter(fem==1) %>% filter(period_new==3 & nr_short!=1), outcome = "Interaction", sample = "REDS - Entire", model = reds.lm2, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm2.pval, coef_dif_p_st = reds.lm2.pval),
  list(data = reds.st.matched %>% filter(fem==1), outcome = "Conflict", sample = "REDS - Entire", model = reds.lm3, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm3.pval, coef_dif_p_st = reds.lm3.pval),
  list(data = reds.scst.matched %>% filter(fem==1), outcome = "Conflict", sample = "REDS - Entire", model = reds.lm4, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm4.pval, coef_dif_p_st = reds.lm4.pval)
  
)

# Create tables and plots of net effects
plots = lapply(effect_list, FUN = function(x){
  effect_size_plot(data = x$data, outcome = x$outcome, sample = x$sample, model = x$model, 
                   vstvar = x$vstvar, vwvar = x$vwvar, coef_dif_p_w = x$coef_dif_p_w,
                   coef_dif_p_st = x$coef_dif_p_st)
})

## Net effect tables
effect_tables = bind_rows(plots[[1]] %>% mutate(outcome="Marriage") %>% mutate(quota = "ST"),
                          plots[[2]] %>% mutate(outcome="Marriage") %>% mutate(quota = "ST/SC"),
                          plots[[3]] %>% mutate(outcome="Interaction") %>% mutate(quota = "ST"),
                          plots[[4]] %>% mutate(outcome="Interaction") %>% mutate(quota = "ST/SC"),
                          plots[[5]] %>% mutate(outcome="Conflict") %>% mutate(quota = "ST"),
                          plots[[6]] %>% mutate(outcome="Conflict") %>% mutate(quota = "ST/SC"))






####### Women only figure ######
effect_tables %>%
  mutate(out_fac = case_when(
    outcome=="Interaction" ~ 1,
    outcome=="Marriage" ~ 2,
    outcome=="Conflict" ~ 3
  )) %>%
  mutate(out_fac = factor(out_fac, labels = c("Interaction", "Marriage", "Conflict"))) %>%
  mutate(lab = as.character(lab_order)) %>%
  filter(lab=="Women's Res."|lab=="Combined Women Effect") %>%
  mutate(lab = ifelse(lab=="Combined Women Effect", "Women's Quota \n When Caste Quota=1", lab)) %>%
  mutate(lab = ifelse(lab=="Women's Res.", "Women's Quota \n When Caste Quota=0", lab)) %>%
  mutate(labels = ifelse(lab_order=="Women's Res.", 2, 1)) %>%
  mutate(lab_fac = factor(labels, labels = c("Two-Dimensional\nWomen's Quota\n(Caste Quota=1)", "One-Dimensional\nWomen's Quota\n(Caste Quota=0)"))) %>%
  ggplot(., aes(x = lab_fac, y = estimate, color = quota))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray48")+
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate + 1.96*std.error), size = 1.5,  position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error), size = 3, position = position_dodge(width = 0.4))+
  geom_pointrange(size = 1.5, aes(shape = quota, ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error),
                  position = position_dodge(width = 0.4), fill = "white") +
  geom_text(aes(label = round(estimate,2)), vjust = 2, size = 6.5,
            position = position_dodge(width = 0.4), show.legend = F)+
  scale_color_manual(values = c("#e66101", "#b2abd2"), name = "Caste Quota Analyzed")+
  scale_shape_manual(values = c(21, 24), name = "Caste Quota Analyzed")+
  coord_flip()+
  facet_grid(.~out_fac, scales = "free")+
  labs(x = "", y = "Estimate (95% CI)")+
  theme_pubr()+
  theme(strip.text = element_text(size = 20),
        axis.text = element_text(size = 20),
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 20))

ggsave(paste0(fig.out, "FigureE8.pdf"),
       width = 12, height = 7)



####### Combined quotas (ST women) ######
# INTER-CASTE MARRIAGE
controls_dn = "age"
base_dn = paste("marry_other_caste_bin ~ vwreserv*vstreserv")
dn.lm1 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st) %>%
              filter(fem==1 & st==1),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata")

dn.lm1.pval = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lh_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st) %>%
              filter(fem==1 & st==1),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code,
            se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv")
dn.lm1.pval = dn.lm1.pval$lh$p.value


# INTER-CASTE MARRIAGE - SC/ST
base_dn = paste("marry_other_caste_bin ~ vwreserv*vscstreserv")
dn.lm2 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst) %>%
              filter(fem==1 & (st==1 | sc==1)),
            fixed_effects = ~uniq_dist_code, clusters = uniq_panch_code, se_type = "stata")

dn.lm2.pval = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lh_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst) %>%
              filter(fem==1 & (st==1 | sc==1)),
            fixed_effects = ~uniq_dist_code, 
            clusters = uniq_panch_code, se_type = "stata",
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")
dn.lm2.pval = dn.lm2.pval$lh$p.value




# EASE OF INTER-CASTE INTERACTIONS
base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
controls = "top20land + birthyr + vwreserv_gp12 + vstreserv_gp12"
reds.lm1 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1) %>%
              filter(fem==1 & Caste_ST==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm1.pval = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1) %>%
              filter(fem==1 & Caste_ST==1),
            fixed_effects = ~districtid, clusters = villageid,
            se_type = "stata", linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv") 
reds.lm1.pval = reds.lm1.pval$lh$p.value

# EASE OF INTER-CASTE INTERACTIONS - SC/ST 
controls = "top20land + birthyr + vwreserv_gp12 + vscstreserv_gp12"
base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"
reds.lm2 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1) %>%
              filter(fem==1 & (Caste_ST==1|Caste_SC==1)),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")

reds.lm2.pval = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lh_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1) %>%
              filter(fem==1 & (Caste_ST==1|Caste_SC==1)),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata",
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")
reds.lm2.pval = reds.lm2.pval$lh$p.value


# # CASTE-BASED CONFLICT
base = "caste_based_inter_family_disputes ~ vwreserv*vstreserv"
controls = "top20land + birthyr + vwreserv_prior + vstreserv_prior + prop_STnow"
reds.lm3 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(fem==1 & Caste_ST==1), fixed_effects = ~period_new + districtid, clusters = villageid, se_type = "stata")

reds.lm3.pval = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.st.matched %>% filter(fem==1 & Caste_ST==1),
            fixed_effects = ~period_new + districtid,
            clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vwreserv = vwreserv + vwreserv:vstreserv")
reds.lm3.pval = reds.lm3.pval$lh$p.value




# # CASTE-BASED CONFLICT
base = "caste_based_inter_family_disputes ~ vwreserv*vscstreserv"
controls = "top20land + birthyr + vwreserv_prior + vscstreserv_prior + prop_STnow + prop_SCnow"
reds.lm4 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(fem==1 & (Caste_ST==1 | Caste_SC==1)), fixed_effects = ~period_new + districtid, clusters = villageid, se_type = "stata")

reds.lm4.pval = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.scst.matched %>% filter(fem==1 & (Caste_ST==1 | Caste_SC==1)),
            fixed_effects = ~period_new + districtid, clusters = villageid,
            se_type = "stata", linear_hypothesis = "vwreserv = vwreserv + vwreserv:vscstreserv")
reds.lm4.pval = reds.lm4.pval$lh$p.value



# Create dataset of models to loop over
effect_list = list(
  list(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st) %>% filter(fem==1 & st==1), outcome = "Marriage", sample = "DN - Entire", model = dn.lm1, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = dn.lm1.pval, coef_dif_p_st = dn.lm1.pval),
  list(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_scst) %>% filter(fem==1 & (st==1|sc==1)), outcome = "Marriage", sample = "DN - Entire", model = dn.lm2, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = dn.lm2.pval, coef_dif_p_st = dn.lm2.pval),
  list(data = reds.st.matched %>% filter(fem==1 & Caste_ST==1) %>% filter(period_new==3 & nr_short!=1), outcome = "Interaction", sample = "REDS - Entire", model = reds.lm1, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm1.pval, coef_dif_p_st = reds.lm1.pval),
  list(data = reds.scst.matched %>% filter(fem==1 & (Caste_ST==1|Caste_SC==1)) %>% filter(period_new==3 & nr_short!=1), outcome = "Interaction", sample = "REDS - Entire", model = reds.lm2, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm2.pval, coef_dif_p_st = reds.lm2.pval),
  list(data = reds.st.matched %>% filter(fem==1 & Caste_ST==1), outcome = "Conflict", sample = "REDS - Entire", model = reds.lm3, vstvar = "vstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm3.pval, coef_dif_p_st = reds.lm3.pval),
  list(data = reds.scst.matched %>% filter(fem==1 & (Caste_ST==1|Caste_SC==1)), outcome = "Conflict", sample = "REDS - Entire", model = reds.lm4, vstvar = "vscstreserv", vwvar = "vwreserv", coef_dif_p_w = reds.lm4.pval, coef_dif_p_st = reds.lm4.pval))

# Create tables and plots of net effects
plots = lapply(effect_list, FUN = function(x){
  effect_size_plot(data = x$data, outcome = x$outcome, sample = x$sample, model = x$model, 
                   vstvar = x$vstvar, vwvar = x$vwvar, coef_dif_p_w = x$coef_dif_p_w, 
                   coef_dif_p_st = x$coef_dif_p_st)
})

## Net effect tables
effect_tables = bind_rows(plots[[1]] %>% mutate(outcome="Marriage") %>% mutate(quota = "ST"),
                          plots[[2]] %>% mutate(outcome="Marriage") %>% mutate(quota = "ST/SC"),
                          plots[[3]] %>% mutate(outcome="Interaction") %>% mutate(quota = "ST"),
                          plots[[4]] %>% mutate(outcome="Interaction") %>% mutate(quota = "ST/SC"),
                          plots[[5]] %>% mutate(outcome="Conflict") %>% mutate(quota = "ST"),
                          plots[[6]] %>% mutate(outcome="Conflict") %>% mutate(quota = "ST/SC"))





####### ST/SC Women only figure ######
effect_tables %>%
  mutate(out_fac = case_when(
    outcome=="Interaction" ~ 1,
    outcome=="Marriage" ~ 2,
    outcome=="Conflict" ~ 3
  )) %>%
  mutate(out_fac = factor(out_fac, labels = c("Interaction", "Marriage", "Conflict"))) %>%
  mutate(lab = as.character(lab_order)) %>%
  filter(lab=="Women's Res."|lab=="Combined Women Effect") %>%
  mutate(lab = ifelse(lab=="Combined Women Effect", "Women's Quota \n When Caste Quota=1", lab)) %>%
  mutate(lab = ifelse(lab=="Women's Res.", "Women's Quota \n When Caste Quota=0", lab)) %>%
  mutate(labels = ifelse(lab_order=="Women's Res.", 2, 1)) %>%
  mutate(lab_fac = factor(labels, labels = c("Two-Dimensional\nWomen's Quota\n(Caste Quota=1)", "One-Dimensional\nWomen's Quota\n(Caste Quota=0)"))) %>%
  ggplot(., aes(x = lab_fac, y = estimate, color = quota))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray48")+
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate + 1.96*std.error), size = 1.5,  position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error), size = 3, position = position_dodge(width = 0.4))+
  geom_pointrange(size = 1.5, aes(shape = quota, ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error),
                  position = position_dodge(width = 0.4), fill = "white") +
  geom_text(aes(label = round(estimate,2)), vjust = 2, size = 6,
            position = position_dodge(width = 0.4), show.legend = F)+
  scale_color_manual(values = c("#e66101", "#b2abd2"), name = "Caste Quota Analyzed")+
  scale_shape_manual(values = c(21, 24), name = "Caste Quota Analyzed")+
  coord_flip()+
  facet_grid(.~out_fac, scales = "free")+
  labs(x = "", y = "Estimate (95% CI)")+
  theme_pubr()+
  theme(strip.text = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.text.x = element_text(size = 16),
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 20))

ggsave(paste0(fig.out, "FigureE9.pdf"),
       width = 12, height = 7)


####  Long-term effect of quotas ####
## Individual plots
effect_size_plot_longterm = function(data, outcome, sample, model, vstvar, vwvar, vw_vst_var, coef_dif_p_w, coef_dif_p_st){
  vstvar = ifelse(is.null(vstvar), "vstreserv", vstvar)
  vwvar = ifelse(is.null(vwvar), "vwreserv", vwvar)
  vw_vst_var = ifelse(is.null(vw_vst_var), paste0(vwvar,":",vstvar), vw_vst_var)
  # Calculate effect sizes
  sample = paste0(sample, ", Outcome: ", outcome)
  women_effect = tidy(model)$estimate[tidy(model)$term==vwvar]
  st_effect = tidy(model)$estimate[tidy(model)$term==vstvar]
  form_st = paste(vstvar, " + ", vw_vst_var, " = 0")
  combined_st_lht = glht(model, linfct = form_st)
  combined_st_effect = tidy(combined_st_lht)$estimate
  
  form_w = paste(vwvar, " + ", vw_vst_var, " = 0")
  combined_w_lht = glht(model, linfct = form_w)
  combined_w_effect = tidy(combined_w_lht)$estimate
  
  
  # Add standard errors
  women_se = tidy(model)$std.error[tidy(model)$term==vwvar]
  st_se = tidy(model)$std.error[tidy(model)$term==vstvar]
  combined_w_se = tidy(combined_w_lht)$std.error
  combined_st_se = tidy(combined_st_lht)$std.error
  
  # Add p-values
  women_p = tidy(model)$p.value[tidy(model)$term==vwvar]
  st_p = tidy(model)$p.value[tidy(model)$term==vstvar]
  combined_w_p = tidy(combined_w_lht)$adj.p.value
  combined_st_p = tidy(combined_st_lht)$adj.p.value
  
  
  figure_table = data.frame(
    estimate = c(women_effect, st_effect, combined_w_effect, combined_st_effect),
    std.error = c(women_se, st_se, combined_w_se, combined_st_se),
    p.value = c(women_p, st_p, combined_w_p, combined_st_p),
    coef_dif_p = c(coef_dif_p_w, coef_dif_p_st,NA_real_, NA_real_),
    labels = c(4, 3, 2, 1)
  ) %>% mutate(lab_order = factor(labels, labels = c("Combined ST Effect", "Combined Women Effect", "ST Res.","Women's Res.")))
  

  return(figure_table)
}
base = "caste_interact_bi_comp ~ vwreserv_gp12 + vstreserv_gp12 + vwrXvst_gp12"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + prop_STnow"
reds.lm5 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

reds.lm5.pval =  paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata", 
            linear_hypothesis = "vwreserv_gp12 = vwreserv_gp12 + vwrXvst_gp12") 
reds.lm5.pval = reds.lm5.pval$lh$p.value

base = "caste_interact_bi_comp ~ vwreserv_gp12 + vscstreserv_gp12 + vwrXvscst_gp12 + prop_STnow + prop_SCnow"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC"
reds.lm7 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")

reds.lm7.pval = paste(base, "+", controls) %>% as.formula() %>%
  lh_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata",
            linear_hypothesis = "vwreserv_gp12 = vwreserv_gp12 + vwrXvscst_gp12")






reds.lm7.pval = reds.lm7.pval$lh$p.value

# Create dataset of models to loop over
effect_list = list(
  list(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1), outcome = "Interaction", sample = "REDS - Entire", 
       model = reds.lm5, vstvar = "vstreserv_gp12", vwvar = "vwreserv_gp12", vw_vst_var = "vwrXvst_gp12",
       coef_dif_p_w = reds.lm5.pval, coef_dif_p_st = reds.lm5.pval),
  list(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1), outcome = "Interaction", sample = "REDS - Entire", 
       model = reds.lm7, vstvar = "vscstreserv_gp12", vwvar = "vwreserv_gp12", vw_vst_var = "vwrXvscst_gp12",
       coef_dif_p_w = reds.lm7.pval, coef_dif_p_st = reds.lm7.pval))

# Create tables and plots of net effects
plots = lapply(effect_list, FUN = function(x){
  effect_size_plot_longterm(data = x$data, outcome = x$outcome, sample = x$sample, model = x$model, 
                   vstvar = x$vstvar, vwvar = x$vwvar, vw_vst_var = x$vw_vst_var,
                   coef_dif_p_w = x$coef_dif_p_w, coef_dif_p_st = x$coef_dif_p_st)
})

## Net effect tables
effect_tables = bind_rows(plots[[1]] %>% mutate(outcome="Interaction") %>% mutate(quota = "ST"),
                          plots[[2]] %>% mutate(outcome="Interaction") %>% mutate(quota = "ST/SC"))
effect_tables = effect_tables %>% 
  mutate(coef_dif_p_chr = ifelse(!is.na(coef_dif_p), paste0("Coef. test p-val=0.000"), NA_character_))





####### Women only figure ######
effect_tables %>% 
  filter(outcome=="Interaction") %>%
  mutate(lab = as.character(lab_order)) %>%
  filter(lab=="Women's Res."|lab=="Combined Women Effect") %>%
  mutate(lab = ifelse(lab=="Combined Women Effect", "Women's Quota \n When Caste Quota=1", lab)) %>%
  mutate(lab = ifelse(lab=="Women's Res.", "Women's Quota \n When Caste Quota=0", lab)) %>%
  mutate(labels = ifelse(lab_order=="Women's Res.", 2, 1)) %>%
  mutate(lab_fac = factor(labels, labels = c("Prior Two-Dimensional\nWomen's Quota\n(Caste Quota=1)", "Prior One-Dimensional\nWomen's Quota\n(Caste Quota=0)"))) %>%
  ggplot(., aes(x = lab_fac, color = quota))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray48")+
  geom_linerange(aes(ymin = estimate-1.96*std.error, ymax = estimate + 1.96*std.error), size = 1.5,  position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error), size = 3, position = position_dodge(width = 0.4))+
  geom_pointrange(size = 1.5, aes(y = estimate, shape = quota, ymin = estimate-1.64*std.error, ymax = estimate + 1.64*std.error),
                  position = position_dodge(width = 0.4), fill = "white") + 
  geom_text(aes(label = round(estimate,2), y = estimate,), vjust = 2, size = 6.5,
            position = position_dodge(width = 0.4), show.legend = F)+
  geom_text(aes(y = 0.006, label = coef_dif_p_chr), size = 5.5, hjust = 0,
            position = position_dodge(width = 0.4), show.legend = F)+
  scale_color_manual(values = c("#e66101", "#b2abd2"), name = "Caste Quota Analyzed")+
  scale_shape_manual(values = c(21, 24), name = "Caste Quota Analyzed")+
  coord_flip()+
  facet_grid(.~outcome, scales = "free")+
  labs(x = "", y = "Estimate (95% CI)")+
  theme_pubr()+
  theme(strip.text = element_text(size = 22), 
        axis.text = element_text(size = 22), 
        axis.title = element_text(size = 22),
        legend.text = element_text(size = 22),
        legend.title = element_text(size = 22),
        axis.title.y =element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())
ggsave(paste0(fig.out, "Figure1b.pdf"),
       width = 10, height = 6.5)



